#!/usr/bin/perl
#===============================================================================
# foldTimeSeries.pl: script to fold a time series
#===============================================================================
# Last update: 31 July 2008
#===============================================================================


#use warnings;
#use strict;

use Math::BigFloat;
use POSIX;

# command line options
print STDERR "============================================\n";
print STDERR "             foldTimeSeries.pl\n";
print STDERR "============================================\n";
$argc = $#ARGV+1;
if ($argc != 4) {
    print STDERR "ERROR: incorrect number of command line args ($argc/4)\n";
    print STDERR "ERROR:\n";
    print STDERR "USAGE: foldTimeSeries.pl [data file] [fold period] [n bins] [output file root name]\n";
    print STDERR "\n";
    exit(99);
}

# save command line arguments and initilise/declare variables
$dataFile   = $ARGV[0]; 
$foldPeriod = Math::BigFloat->new($ARGV[1]); 
$timeOffset = 0.0;
$outputFile = $ARGV[3];
$nBins      = $ARGV[2];                      # number of phase bins
$binWidth   = 1.0/$nBins;               # bin width in phase units (0->2pi)
$nData      = 0;                        # number of data points
$phase      = Math::BigFloat->new(0.0); # phase of the signal (in units of 2pi)
$foldFreq   = 1.0/$foldPeriod;
$epoch      = Math::BigFloat->new($timeOffset);   # time origin offset
#$epoch     = 2450023.62 - 485.*$p;
$outputFile = "${outputFile}-p${foldFreq}-n${nBins}.dat";

# read the data file
$i = 0;
open(DATA,"$dataFile") || die "cannot open: $dataFile\n\n";
while(<DATA>) {
    split;
    $time[$i] = Math::BigFloat->new($_[0]);
    $amp[$i] = Math::BigFloat->new($_[1]);
    $a = Math::BigFloat->new(0.0);
    $b = Math::BigFloat->new(0.0);
    $a =  ($time[$i]-$epoch)/$foldPeriod;
    $b =  floor(($time[$i]-$epoch)/$foldPeriod);
    $phase[$i] = ($a - $b);
    printf(STDERR "(%3d)  time: %.12f\tamp: %.12f\tphase: %.18f  ($phase[$i])\n",$i,$time[$i],$amp[$i],$phase[$i]);
    $i++;
}
close(DATA);
$nData = $i;
 

# print a summary of the input options
print STDERR "--------------- INPUT SUMMARY --------------\n";
print STDERR "Datafile        = $dataFile\n";
print STDERR "nData           = $nData\n";
print STDERR "\n";
print STDERR "Folding freq    = $foldFreq\n";
print STDERR "Folding period  = $foldPeriod\n";
print STDERR "Time offset     = $timeOffset\n";
print STDERR "\n";
print STDERR "# bins          = $nBins\n";
print STDERR "Bin width       = $binWidth (2pi)\n";
print STDERR "\n";
print STDERR "Output file     = $outputFile\n";
print STDERR "--------------------------------------------\n";

# put data into bins
$totalBinned = 0;
$binOffset = Math::BigFloat->new($binWidth/2.0);
print STDERR "--------------------------------------------\n";
for ($iBin = 0;$iBin < $nBins; $iBin++) {
    $binStart[$iBin] = Math::BigFloat->new(0.0);
    $binCentre[$iBin]= Math::BigFloat->new(0.0);
    $binEnd[$iBin]   = Math::BigFloat->new(0.0);
    $binStart[$iBin] = $iBin * $binWidth;
    $binCentre[$iBin]= $binStart[$iBin] + $binOffset;
    $binEnd[$iBin]   = $binStart[$iBin] + $binWidth;
    $binCount[$iBin] = 0;
    $binTotal[$iBin] = Math::BigFloat->new(0.0);
    $binMean[$iBin]  = Math::BigFloat->new(0.0);
    printf(STDERR "[bin %2d]\t($binStart[$iBin] <= x < $binEnd[$iBin])\n",$iBin);
    for ($i=0;$i<$nData;$i++) {
        if ( ($phase[$i] >= $binStart[$iBin]) && ($phase[$i] < $binEnd[$iBin]) ) {
            $binCount[$iBin] += 1;
            $binTotal[$iBin] += $amp[$i];
            $totalBinned += 1;
            printf(STDERR "   ($i) phase = $phase[$i]\t\t(amp = %.3f total = %.3f)\n",$amp[$i],$binTotal[$iBin]);
        }
    }
    $binMean[$iBin] = $binTotal[$iBin]/$binCount[$iBin];
    printf STDERR "[bin %2d]\t(Found $binCount[$iBin] points)\n\n",$iBin; 
}
if ($totalBinned != $nData) {
    print STDERR "ERROR: not all the data points were binned!\n";
    exit(1);
}

# caluclate mean amplitude and standard deviation in each bin
for ($iBin=0; $iBin<$nBins;$iBin++) { 
    $binStdDev[$iBin] = Math::BigFloat->new(0.0);
    for ($i=0;$i<$nData;$i++) {
        if ( ($phase[$i] >= $binStart[$iBin]) && ($phase[$i] < $binEnd[$iBin]) ) {
            $sumSqDev[$iBin] += ($amp[$i] - $binMean[$iBin])**2;
        }
    }
    $binStdDev[$iBin] = sqrt($sumSqDev[$iBin]/$binCount[$iBin]);
    printf(STDERR "($iBin) count = $binCount[$iBin] mean = %.5f sumSq = %.5f stddev = %.5f\n",$binMean[$iBin],$sumSqDev[$iBin],$binStdDev[$iBin]); 
}

open(OUT,">$outputFile") || die "cannot open: $outputFile\n\n";
print STDERR "--------------------------------------------\n";
#========== print output to STDOUT ===========================   
print OUT "# phase | mean | stdDev | #points | iBin \n";
print STDERR "# phase | mean | stdDev | #points | iBin \n";
for ($iBin=0;$iBin<$nBins;$iBin++) { 
    printf(STDOUT "%f %f %f %d %d\n",$binCentre[$iBin],$binMean[$iBin],$binStdDev[$iBin],$binCount[$iBin],$iBin);
    printf(OUT "%f %f %f %d %d\n",$binCentre[$iBin],$binMean[$iBin],$binStdDev[$iBin],$binCount[$iBin],$iBin);
    
}    
print STDERR "--------------------------------------------\n";

close(OUT);
